%******************************************************************
%** This program is specifically written to replicate Model I in **
%** Hellwig and Samuels's (2007) CPS article.                    **
%** This program applies the revised TRM model with              **
%** constrained optimization. (TRMCO)                            **
%** Created by Min-Hua Huang, Dec 19, 2011                       **
%******************************************************************
clear; clc; format long;
load Hellwig1.mat; load stataresult1.mat; 
% data manipulaitons
a=0; b=100; y=Hellwig1(:,1); n=size(y,1); xr=Hellwig1(:,2:12);
mf=size(xr,2); minx=min(xr)';
for i=1:1:mf
    xr(:,i)=xr(:,i)-minx(i);
end
cr1=xr(:,2).*xr(:,3); cr2=xr(:,2).*xr(:,4);
x=[xr(:,1:3) cr1 xr(:,4) cr2 xr(:,5:11)];
code=Hellwig1(:,13);
% save the data for identifying boundary violations (identify.m) 
save x.mat x;
% parameter setups for TRMCO
m=size(x,2); dmy=4; runmax=101; tau=14; tolv=10^-4;
trmcom=trmcod(y,x,[trmb;trmsig],a,b,tau,runmax,dmy,tolv);
%**** TRMCO ****
trmcob=trmcom(1,1:m+1)';
trmcose=robuster(y,[ones(n,1) x],trmcom(1,1:m+2)',code);
trmcoyhat=[ones(n,1) x]*trmcob;
trmcoymax=max(trmcoyhat); 
trmcoymin=min(trmcoyhat);
trmcot=trmcob./trmcose(1:m+1);
LLtrmco=trmcom(m+5,2);
viotemp1=identify(trmcom(1,1:m+2)');
trmcovios=(2*m+6)*2-sum(sum(viotemp1<(b+10^-4)))-sum(sum(viotemp1>a-10^-4));
%**** OLS ****
olsyhat=[ones(n,1) x]*olsb;
olsLL=LL(y,[ones(n,1) x],[olsb;olssig],a,b);
max1=max(x)';
min1=min(x)';
xlimit=[[1;max1] [1;min1]];
v=vsign(olsb);
olsci=cpar([olsb;olssig],a,b,[ones(n,1) x],xlimit,v,dmy);
viotemp2=identify([olsb;olssig]);
olsvios=(2*m+6)*2-sum(sum(viotemp2<(b+10^-4)))-sum(sum(viotemp2>a-10^-4));
olst=olsb./olsse;
%**** TRM ****
trmyhat=[ones(n,1) x]*trmb;
trmymax=max(trmyhat);
trmymin=min(trmyhat);
LLtrm=LL(y,[ones(n,1) x],[trmb;trmsig],a,b);
v=vsign(trmb);
trmci=cpar([trmb;trmsig],a,b,[ones(n,1) x],xlimit,v,dmy);
viotemp3=identify([trmb;trmsig]);
trmvios=(2*m+6)*2-sum(sum(viotemp3<(b+10^-4)))-sum(sum(viotemp3>a-10^-4));
trmt=trmb./trmse;
% ****  OUTPUT ****
fprintf('    \nReplification of Hellwig and Samuels (2007), Model I \n'  );
fprintf('                          OLS                              \n'  );   
fprintf('                                Coeff        SE       t    \n');   
fprintf('Previous vote                   %6.3f     %6.3f   %6.3f\n',olsb(2),olsse(2),olst(2));
fprintf('Economy                         %6.3f     %6.3f   %6.3f\n',olsb(3),olsse(3),olst(3));
fprintf('Trade Openness                  %6.3f     %6.3f   %6.3f\n',olsb(4),olsse(4),olst(4));
fprintf('Economy*Trade Openness          %6.3f     %6.3f   %6.3f\n',olsb(5),olsse(5),olst(5));
fprintf('Presidential Election           %6.3f     %6.3f   %6.3f\n',olsb(6),olsse(6),olst(6));
fprintf('Economy*Presidential Election   %6.3f     %6.3f   %6.3f\n',olsb(7),olsse(7),olst(7));
fprintf('Re-election                     %6.3f     %6.3f   %6.3f\n',olsb(8),olsse(8),olst(8));
fprintf('Effective Number of Parties     %6.3f     %6.3f   %6.3f\n',olsb(9),olsse(9),olst(9));
fprintf('Income                          %6.3f     %6.3f   %6.3f\n',olsb(10),olsse(10),olst(10));
fprintf('Africa                          %6.3f     %6.3f   %6.3f\n',olsb(11),olsse(11),olst(11));
fprintf('Asia                            %6.3f     %6.3f   %6.3f\n',olsb(12),olsse(12),olst(12));
fprintf('Central and Eastern Europe      %6.3f     %6.3f   %6.3f\n',olsb(13),olsse(13),olst(13));
fprintf('Latin America and the Caribbean %6.3f     %6.3f   %6.3f\n',olsb(14),olsse(14),olst(14));
fprintf('Constant                        %6.3f     %6.3f   %6.3f\n',olsb(1),olsse(1),olst(1));
fprintf('Sigma_hat                       %6.3f                  \n',olssig);
fprintf('Log pseudolikelihood=           %6.3f\n',olsLL);
fprintf('Max(y-hat)=                     %6.3f\n',max(olsyhat));
fprintf('Min(y-hat)=                     %6.3f\n',min(olsyhat));
fprintf('Bourndary Violations=           %6.0f\n',olsvios);
fprintf('N=                              %6.0f\n',n);
% **********************************************************************
fprintf('                          TRM                              \n'     );   
fprintf('                                Coeff        SE       t    \n');   
fprintf('Previous vote                   %6.3f     %6.3f   %6.3f\n',trmb(2),trmse(2),trmt(2));
fprintf('Economy                         %6.3f     %6.3f   %6.3f\n',trmb(3),trmse(3),trmt(3));
fprintf('Trade Openness                  %6.3f     %6.3f   %6.3f\n',trmb(4),trmse(4),trmt(4));
fprintf('Economy*Trade Openness          %6.3f     %6.3f   %6.3f\n',trmb(5),trmse(5),trmt(5));
fprintf('Presidential Election           %6.3f     %6.3f   %6.3f\n',trmb(6),trmse(6),trmt(6));
fprintf('Economy*Presidential Election   %6.3f     %6.3f   %6.3f\n',trmb(7),trmse(7),trmt(7));
fprintf('Re-election                     %6.3f     %6.3f   %6.3f\n',trmb(8),trmse(8),trmt(8));
fprintf('Effective Number of Parties     %6.3f     %6.3f   %6.3f\n',trmb(9),trmse(9),trmt(9));
fprintf('Income                          %6.3f     %6.3f   %6.3f\n',trmb(10),trmse(10),trmt(10));
fprintf('Africa                          %6.3f     %6.3f   %6.3f\n',trmb(11),trmse(11),trmt(11));
fprintf('Asia                            %6.3f     %6.3f   %6.3f\n',trmb(12),trmse(12),trmt(12));
fprintf('Central and Eastern Europe      %6.3f     %6.3f   %6.3f\n',trmb(13),trmse(13),trmt(13));
fprintf('Latin America and the Caribbean %6.3f     %6.3f   %6.3f\n',trmb(14),trmse(14),trmt(14));
fprintf('Constant                        %6.3f     %6.3f   %6.3f\n',trmb(1),trmse(1),trmt(1));
fprintf('Sigma_hat                       %6.3f                  \n',trmsig); 
fprintf('Log pseudolikelihood=           %6.3f\n',LLtrm);
fprintf('Max(y-hat)=                     %6.3f\n',trmymax);
fprintf('Min(y-hat)=                     %6.3f\n',trmymin);
fprintf('Bourndary Violations=           %6.0f\n',trmvios);
fprintf('N=                              %6.0f\n',n);
% **********************************************************************
fprintf('                          TRMCO                            \n'     );   
fprintf('                                Coeff        SE       t    \n');   
fprintf('Previous vote                   %6.3f     %6.3f   %6.3f\n',trmcob(2),trmcose(2),trmcot(2));
fprintf('Economy                         %6.3f     %6.3f   %6.3f\n',trmcob(3),trmcose(3),trmcot(3));
fprintf('Trade Openness                  %6.3f     %6.3f   %6.3f\n',trmcob(4),trmcose(4),trmcot(4));
fprintf('Economy*Trade Openness          %6.3f     %6.3f   %6.3f\n',trmcob(5),trmcose(5),trmcot(5));
fprintf('Presidential Election           %6.3f     %6.3f   %6.3f\n',trmcob(6),trmcose(6),trmcot(6));
fprintf('Economy*Presidential Election   %6.3f     %6.3f   %6.3f\n',trmcob(7),trmcose(7),trmcot(7));
fprintf('Re-election                     %6.3f     %6.3f   %6.3f\n',trmcob(8),trmcose(8),trmcot(8));
fprintf('Effective Number of Parties     %6.3f     %6.3f   %6.3f\n',trmcob(9),trmcose(9),trmcot(9));
fprintf('Income                          %6.3f     %6.3f   %6.3f\n',trmcob(10),trmcose(10),trmcot(10));
fprintf('Africa                          %6.3f     %6.3f   %6.3f\n',trmcob(11),trmcose(11),trmcot(11));
fprintf('Asia                            %6.3f     %6.3f   %6.3f\n',trmcob(12),trmcose(12),trmcot(12));
fprintf('Central and Eastern Europe      %6.3f     %6.3f   %6.3f\n',trmcob(13),trmcose(13),trmcot(13));
fprintf('Latin America and the Caribbean %6.3f     %6.3f   %6.3f\n',trmcob(14),trmcose(14),trmcot(14));
fprintf('Constant                        %6.3f     %6.3f   %6.3f\n',trmcob(1),trmcose(1),trmcot(1));
fprintf('Sigma_hat                       %6.3f                  \n',trmcom(1,m+2)); 
fprintf('Log pseudolikelihood=           %6.3f\n',LLtrmco);
fprintf('Max(y-hat)=                     %6.3f\n',trmcoymax);
fprintf('Min(y-hat)=                     %6.3f\n',trmcoymin);
fprintf('Bourndary Violations=           %6.0f\n',trmcovios);
fprintf('N=                              %6.0f\n',n);
save overall_result1.mat;